if(!require("mlr3verse")) install.packages("mlr3verse")
if(!require("tidyverse")) install.packages("tidyverse")
if(!require("ggplot2")) install.packages("ggplot2")
if(!require("rms")) install.packages("rms")

rm(list = ls())

df_select <- readRDS("~/analysis/lyz_ml/rds/step_03_mlr3_select.RDS")
colnames(df_select) <- c("group",'Bronchoscopy', 'Age', 'Days of illness prior to admission', 'WBC before treatment',
                         'PLR before treatment',"CRP before treatment", 'LDH before treatment') # 
colnames(df_select) <- make.names(names(df_select) )

df_select$group <- as.numeric(df_select$group)-1

# 环境
ddist <- datadist(df_select)
options(datadist='ddist')

logis <- lrm(group ~ ., data=df_select )
nom1 <- nomogram(logis, fun=plogis,
                 fun.at=c(0.001,0.1,0.25,0.5,0.75,0.9,0.99),
                 lp=T, # 是否显示线性概率
                 funlabel="Risk of MP")  

pdf(file = 'result/step_06_logis_nom.pdf',width = 12,height = 6)
plot(nom1) 
dev.off()

# 模型参数
logis <- glm(group ~ ., data = df_select, family = binomial)
summary(logis)$coef
coef(logis)
anova(object = logis, test = "Chisq")
coef(logis)
exp(coef(logis)) # OR


# 模型比较
library("MASS")
step.model <- step(object = logis, trace = 0)
summary(step.model)
anova(object = step.model, test = "Chisq")
anova(logis, step.model, test = "Chisq")


# 模型图
# Fill predicted values using regression model
df_select$pred = predict( logis, df_select , type="response")

ggplot(df_select, aes(x=pred, y=group)) + 
  geom_point(alpha=0.2,size=0.5) +
  stat_smooth(method="glm", color="grey60", se=FALSE,
              method.args = list(family=binomial))+
  ylab("MP")+
  xlab("Predict")+
  geom_point(aes(x = 0.5,y = 0.5),color='red')+
  theme_classic()+
  scale_y_continuous( breaks = c(0,1))+
  scale_x_continuous( breaks = seq(0.1,1,0.2),limits = c(0.1,0.9))

ggsave(filename = 'result/step_06_logistic_s.pdf',width = 14,height = 8,dpi = 300,units = 'cm')


